Script 8 for Kitchel et al.Ā 2024 in prep taxonomic diversity manuscript.
library(data.table)
library(MuMIn)
library(ggplot2)
library(cowplot)
library(lme4)
library(stringr)
library(nlme)
#pull in function to calculate model estimates and standard errors
source(here::here("analysis_code","extract_coefficients_function.R"))
###Predicts annual dissimilarity with annual characteristics, temperature and fishing values
Pull in - region areas (if not already loaded) - region characteristics (if not already loaded); saveRDS(FishGlob_richness_year_survey, file = here::here(āoutputā,āFishGlob_richness_year_survey.Rdsā)) - fishing (if not already loaded) - temp (if not already loaded)
Add survey area to dissimilarities data table
#physical area by year
region_area_byyear <- fread(here::here("output","region_area_byyear.csv"))
#merged fishing, temp, dissimilarities
dissimilarities_temp_fishing <- fread(here::here("output","dissimilarities_temp_fishing.csv"))
#combine
dissimilarities_temp_fishing_area <- dissimilarities_temp_fishing[region_area_byyear, on = c("survey_unit","year")]
#only occurrence-based jaccard for these analyses
dissimilarities_temp_fishing_area.jaccard <- dissimilarities_temp_fishing_area[dissimilarity_metric == "jaccard_dissimilarity_index_binary",]
Add in season/Julian day to dissimilarities data table
#load up julian days
dates_regions <- readRDS(here::here("output","dates_regions.rds"))
#most common season per year
# Function to get the most frequent value
get_mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
# Apply the function to each group
dates_regions.r <- dates_regions[, .(season = get_mode(season)), by = c("survey_unit","year")]
#merge
dissimilarities_temp_fishing_area.jaccard <- dates_regions.r[dissimilarities_temp_fishing_area.jaccard, on = c("survey_unit","year")]
Pull in palette and name helper
source(here::here("analysis_code","color_links.R"))
Pull in observed dissimilarity trend values
jaccard_total_coefs.r <- fread(here::here("output","jaccard_total_coefs.r.csv"))
Plot fishing and temperature vs.Ā time for all regions
#######TEMPERATURE
#set order by survey unit for plotting
all_surveys <- levels(as.factor(dissimilarities_temp_fishing_area.jaccard$survey_unit))
setorder(dissimilarities_temp_fishing_area.jaccard, survey_unit)
dissimilarities_temp_fishing_area.jaccard[,Survey_Name_Season:=factor(Survey_Name_Season, levels = unique(dissimilarities_temp_fishing_area.jaccard$Survey_Name_Season), ordered = T)]
(sbt_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
labs(y = "Mean bottom temperature (ĖC)", x = "Year") +
geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
theme_classic() +
theme(axis.text.x = element_text(size = 7)) +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))
ggsave(sbt_time_survey_facet_1_20, path = here::here("figures"), filename = "sbt_time_survey_facet_1_20.jpg", height = 12, width =9)
(sbt_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
labs(y = "Mean bottom temperature (ĖC)", x = "Year") +
geom_point(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), alpha = 0.3) +
geom_smooth(aes(y = as.numeric(yearly_mean_bypoint_avg), x = year), method = "lm") +
scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
theme_classic() +
theme(axis.text.x = element_text(size = 7)) +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))
ggsave(sbt_time_survey_facet_21_34, path = here::here("figures"), filename = "sbt_time_survey_facet_21_34.jpg", height = 12, width =9)
#######FISHING
dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]
(fishing_time_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
labs(y = "Relative fishing catch", x = "Year") +
geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
theme_classic() +
theme(axis.text.x = element_text(size = 7)) +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))
ggsave(fishing_time_survey_facet_1_20, path = here::here("figures"), filename = "fishing_time_survey_facet_1_20.jpg", height = 12, width =9)
(fishing_time_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
labs(y = "Relative fishing catch", x = "Year") +
geom_point(aes(y = summed_tonnes_scaled_byreg, x = year), alpha = 0.3) +
geom_smooth(aes(y = summed_tonnes_scaled_byreg, x = year), method = "lm") +
scale_x_continuous(breaks = ~ axisTicks(., log = FALSE)) +
theme_classic() +
theme(axis.text.x = element_text(size = 7)) +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4))
ggsave(fishing_time_survey_facet_21_34, path = here::here("figures"), filename = "fishing_time_survey_facet_21_34.jpg", height = 12, width =9)
NA
NA
Plot fishing and temperature vs.Ā dissimilarity for all regions
#####MEAN TEMP
(preds_sbt_mean_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
labs(x = "Mean bottom temperature (ĖC)", y = "β-diversity") +
geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_mean_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_1_20.jpg", height = 12, width =9)
(preds_sbt_mean_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
labs(x = "Mean bottom temperature (ĖC)", y = "β-diversity") +
geom_point(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = as.numeric(yearly_mean_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_mean_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_temp_survey_facet_21_34.jpg", height = 12, width =9)
#####MINIMUM TEMP
(preds_sbt_min_temp_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[1:20] & year > 1979]) +
labs(x = "Minimum bottom temperature (ĖC)", y = "β-diversity") +
geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_min_temp_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_1_20.jpg", height = 12, width =9)
(preds_sbt_min_temp_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard[survey_unit %in% all_surveys[21:34] & year > 1979]) +
labs(x = "Minimum bottom temperature (ĖC)", y = "β-diversity") +
geom_point(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = as.numeric(yearly_min_bypoint_avg), y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_min_temp_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_min_temp_survey_facet_21_34.jpg", height = 12, width =9)
#######FISHING
dissimilarities_temp_fishing_area.jaccard.cc <- dissimilarities_temp_fishing_area.jaccard[complete.cases(dissimilarities_temp_fishing_area.jaccard[,summed_tonnes_scaled_byreg]),]
(preds_sbt_mean_fishing_survey_facet_1_20 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[1:20] & year > 1979]) +
labs(x = "Relative fishing catch", y = "β-diversity") +
geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_mean_fishing_survey_facet_1_20, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_1_20.jpg", height = 12, width =9)
(preds_sbt_mean_fishing_survey_facet_21_34 <- ggplot(data = dissimilarities_temp_fishing_area.jaccard.cc[survey_unit %in% all_surveys[c(21:34)] & year > 1979]) +
labs(x = "Relative fishing catch", y = "β-diversity") +
geom_point(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), alpha = 0.3) +
geom_smooth(aes(x = summed_tonnes_scaled_byreg, y = annual_dissimilarity_value), method = "lm") +
facet_wrap(~Survey_Name_Season, scales= "free", ncol = 4) +
theme_classic())
ggsave(preds_sbt_mean_fishing_survey_facet_21_34, path = here::here("figures"), filename = "preds_sbt_mean_fishing_survey_facet_21_34.jpg", height = 12, width =9)
NA
NA
###Plot number of tows per year per region
tows_year <- unique(dissimilarities_temp_fishing_area.jaccard[,.(survey_unit, haul_id_count_annual, area_km, year, Survey_Name_Season)])
min_haul_density <- min(tows_year$haul_id_count_annual/tows_year$area_km)
max_haul_density <- max(tows_year$haul_id_count_annual/tows_year$area_km)
(tow_density_survey_facet_1_20 <- ggplot(data = tows_year[survey_unit %in% all_surveys[1:20]]) +
labs(x = "Year", y = expression("Tow density (tows per km"^2*")")) +
geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
facet_wrap(~Survey_Name_Season, ncol = 4) +
theme_classic())
ggsave(tow_density_survey_facet_1_20, path = here::here("figures"), filename = "tow_density_survey_facet_1_20.jpg", height = 12, width =9)
(tow_density_survey_facet_21_34 <- ggplot(data = tows_year[survey_unit %in% all_surveys[21:34]]) +
labs(x = "Year", y = expression("Tow density (tows per km"^2*")")) +
geom_point(aes(x = year, y = haul_id_count_annual/area_km)) +
ylim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
facet_wrap(~Survey_Name_Season, ncol = 4) +
theme_classic())
ggsave(tow_density_survey_facet_21_34, path = here::here("figures"), filename = "tow_density_survey_facet_21_34.jpg", height = 12, width =9)
#minimum and maximum dissimilarity value
min_dissimilarity <- min(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)
max_dissimilarity <- max(dissimilarities_temp_fishing_area.jaccard$annual_dissimilarity_value)
#how does tow density vary with dissimilarity?
tow_density_dissimilarity_1_20 <-
ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary" & survey_unit %in% all_surveys[1:20]]) +
geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
viridis::scale_color_viridis() +
facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
# xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
# ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
ggsave(tow_density_dissimilarity_1_20, path = here::here("figures"), filename = "tow_density_dissimilarity_1_20.jpg", height = 12, width =9)
tow_density_dissimilarity_21_34 <-
ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary" & survey_unit %in% all_surveys[21:34]]) +
geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = year), size = 2) +
viridis::scale_color_viridis() +
facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
# xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
# ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
ggsave(tow_density_dissimilarity_21_34, path = here::here("figures"), filename = "tow_density_dissimilarity_21_34.jpg", height = 12, width =9)
#Plot all at once?
tow_density_dissimilarity <-
ggplot(data = dissimilarities_temp_fishing_area.jaccard[dissimilarity_metric == "jaccard_dissimilarity_index_binary"]) +
geom_point(aes(x = (haul_id_count_annual/area_km), y = annual_dissimilarity_value, color = Survey_Name_Season), size = 2) +
scale_color_manual(values = color_link$hex) +
labs(x = expression("Tow density (tows per km"^2*")"), y = "β-diversity", color = "Year") +
# xlim(c(min_haul_density-0.0001, max_haul_density+0.0001)) +
# ylim(c(min_dissimilarity-0.01,max_dissimilarity+0.01)) +
theme_classic()
#and by numbers?
result <- dissimilarities_temp_fishing_area.jaccard[, .(max_density = max(haul_id_count_annual/area_km, na.rm = TRUE),
min_density = min(haul_id_count_annual/area_km, na.rm = TRUE)),
by = Survey_Name_Season]
result[,range:= max_density-min_density]
ggplot(data = result) +
geom_point(aes(x = Survey_Name_Season, y = range)) +
theme_classic() +
labs(x = "",y = "Range of tow density values in tow/km^2") +
theme(axis.text.x = element_text(angle = 45, hjust=1))
Plot dissimilarity by season of sampling
dissimilarities_temp_fishing_area.jaccard[,season := factor(season, levels = c("Spring","Summer","Fall","Winter"))]
ggplot(dissimilarities_temp_fishing_area.jaccard) +
geom_boxplot(aes(x = season, y = annual_dissimilarity_value)) +
theme_classic()
###Set up dredge to identify best performing models
First, make data table of model covariates
options(na.action = "na.fail")
dissimilarity_covariates_dredge.dt <- dissimilarities_temp_fishing_area.jaccard[,.
(year, survey_unit,
yearly_mean_bypoint_avg, yearly_max_bypoint_avg, yearly_min_bypoint_avg,yearly_seas_bypoint_avg,
yearly_mean_bypoint_SD, yearly_max_bypoint_SD, yearly_min_bypoint_SD,yearly_seas_bypoint_SD,
yearly_mean_bypoint_avg.s, yearly_max_bypoint_avg.s, yearly_min_bypoint_avg.s,yearly_seas_bypoint_avg.s,
annual_dissimilarity_value,
haul_id_count_annual,
spp_count_annual, depth_annual_avg,
depth_annual_range, latitude_annual_avg,
latitude_annual_range, area_km, season, summed_tonnes_scaled_byreg)]
#merge in with colors for plotting predictions by survey
dissimilarity_covariates_dredge.dt <- color_link[dissimilarity_covariates_dredge.dt, on = "survey_unit"]
#If NA for any covariate, delete row
#View(dissimilarity_covariates_dredge.dt)
#Deleted:
#Before 1980 and after 2019
#Gulf of Saint Laurence South (no depth data)
#No clear SAU match for Rockall Plateau
dissimilarity_covariates_dredge.dt <- dissimilarity_covariates_dredge.dt[complete.cases(dissimilarity_covariates_dredge.dt)]
#Scale and center variables that are not yet scaled and centered
dissimilarity_covariates_dredge.dt[, yearly_mean_bypoint_avg.scaledacrossall := scale(yearly_mean_bypoint_avg)][, yearly_mean_bypoint_SD.scaledacrossall := scale(yearly_mean_bypoint_SD)][, yearly_min_bypoint_avg.scaledacrossall := scale(yearly_min_bypoint_avg)][, yearly_max_bypoint_avg.scaledacrossall := scale(yearly_max_bypoint_avg)][, yearly_seas_bypoint_avg.scaledacrossall := scale(yearly_seas_bypoint_avg)][,haul_id_count_annual.scaledacrossall := scale(haul_id_count_annual)][,spp_count_annual.scaledacrossall := scale(spp_count_annual)][,depth_annual_avg.scaledacrossall := scale(depth_annual_avg)][,depth_annual_range.scaledacrossall := scale(depth_annual_range)] [,latitude_annual_avg.scaledacrossall := scale(latitude_annual_avg)][,latitude_annual_range.scaledacrossall := scale(latitude_annual_range)][,area_km.scaledacrossall := scale(area_km)]
For temperature, we will compare performance of: -mean (scaled across all regions) -max (scaled across all regions) -min (scaled across all regions) -seas (scaled across all regions) -SD (scaled across all regions)
Comparing temp variables AND the presence/absence of a temporal autocorrelation term
#mean temperature
global_mod_mean_temp_lm <- lm(annual_dissimilarity_value ~
survey_unit*yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
data = dissimilarity_covariates_dredge.dt)
#with temporal autocorrelation
global_mod_mean_temp_gls <- gls(annual_dissimilarity_value ~
survey_unit*yearly_mean_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
correlation = corAR1(form=~year | survey_unit),
data = dissimilarity_covariates_dredge.dt)
#maximum temperature
global_mod_max_temp_lm <- lm(annual_dissimilarity_value ~
survey_unit*yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
data = dissimilarity_covariates_dredge.dt)
#with temporal autocorrelation
global_mod_max_temp_gls <-gls(annual_dissimilarity_value ~
survey_unit*yearly_max_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
correlation = corAR1(form=~year | survey_unit),
data = dissimilarity_covariates_dredge.dt)
#minimum temperature
global_mod_min_temp_lm <- lm(annual_dissimilarity_value ~
survey_unit*yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
data = dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt[, survey_unit := factor(survey_unit, levels = sort(all_surveys))]
#with temporal autocorrelation
global_mod_min_temp_gls <-gls(annual_dissimilarity_value ~
survey_unit*yearly_min_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
correlation = corAR1(form=~year | survey_unit),
data = dissimilarity_covariates_dredge.dt)
#temperature seasonality
global_mod_seas_temp_lm <- lm(annual_dissimilarity_value ~
survey_unit*yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
data = dissimilarity_covariates_dredge.dt)
#with temporal autocorrelation
global_mod_seas_temp_gls <- gls(annual_dissimilarity_value ~
survey_unit*yearly_seas_bypoint_avg.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
correlation = corAR1(form=~year | survey_unit),
data = dissimilarity_covariates_dredge.dt)
#standard deviation
global_mod_SD_temp_lm <- lm(annual_dissimilarity_value ~
survey_unit*yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
data = dissimilarity_covariates_dredge.dt)
#with temporal autocorrelation
global_mod_SD_temp_gls <- gls(annual_dissimilarity_value ~
survey_unit*yearly_mean_bypoint_SD.scaledacrossall + #temp and survey unit (possible interaction)
survey_unit*summed_tonnes_scaled_byreg + #fishing effort and survey unit (possible interaction)
area_km.scaledacrossall + #area
latitude_annual_range.scaledacrossall + #latitude range
latitude_annual_avg.scaledacrossall + #latitude avg
depth_annual_range.scaledacrossall + #depth range
depth_annual_avg.scaledacrossall + #depth avg
spp_count_annual.scaledacrossall + #spp #
haul_id_count_annual.scaledacrossall +
season,
correlation = corAR1(form=~year | survey_unit),
data = dissimilarity_covariates_dredge.dt)
#Compare AICcs
View(AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls))
#build data table to report AICc
global_mod_temp_table <- data.table(
`Temporal autocorrelation` = c(
c(rep(F,5),rep(T,5))
),
`Temperature variable` = rep(c(
"Average mean SBT",
"Average maximum SBT",
"Average minimum SBT",
"Average SBT seasonality",
"SBT SD"),2),
deltaAICc = signif(min(AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls)[,2])-AICc(global_mod_mean_temp_lm, global_mod_max_temp_lm, global_mod_min_temp_lm, global_mod_seas_temp_lm, global_mod_SD_temp_lm,
global_mod_mean_temp_gls, global_mod_max_temp_gls, global_mod_min_temp_gls, global_mod_seas_temp_gls, global_mod_SD_temp_gls)[,2],2))
#order by aicc
setorder(global_mod_temp_table,cols = -"deltaAICc")
global_mod_temp_table[,Rank := seq(1,10,by = 1)]
global_mod_sbt_table <- global_mod_temp_table[,.(Rank,`Temporal autocorrelation`,`Temperature variable`, deltaAICc)]
fwrite(global_mod_sbt_table, here::here("output","global_mod_sbt_table.csv"))
####Best performing global model includes minimum temperature (centered and scaled) (Still the case as of July 17, 2024)
Now, look at different combinations of all predictors (min temp) using dredge
Global model: global_mod_min_temp
options(na.action = "na.fail") # prevent fitting sub-models to different datasets
dd <- dredge(global_mod_min_temp_lm)
Fixed term is "(Intercept)"
dd.dt <- data.table(dd)
View(dd)
#only models less than 2 delta AICc (2 models)
dd.dt.2 <- dd.dt[delta <= 2,]
colnames(dd.dt.2)
[1] "(Intercept)"
[2] "area_km.scaledacrossall"
[3] "depth_annual_avg.scaledacrossall"
[4] "depth_annual_range.scaledacrossall"
[5] "haul_id_count_annual.scaledacrossall"
[6] "latitude_annual_avg.scaledacrossall"
[7] "latitude_annual_range.scaledacrossall"
[8] "season"
[9] "spp_count_annual.scaledacrossall"
[10] "summed_tonnes_scaled_byreg"
[11] "survey_unit"
[12] "yearly_min_bypoint_avg.scaledacrossall"
[13] "summed_tonnes_scaled_byreg:survey_unit"
[14] "survey_unit:yearly_min_bypoint_avg.scaledacrossall"
[15] "df"
[16] "logLik"
[17] "AICc"
[18] "delta"
[19] "weight"
#in this step, we delete coefficient values because we will pull them back in later when we calculate both SE and coefficients (other than interaction, which we keep here)
dd.dt.2.formatted <- dd.dt.2[,Rank := as.numeric(rownames(dd.dt.2))][,.(Rank,
`summed_tonnes_scaled_byreg:survey_unit`,
`survey_unit:yearly_min_bypoint_avg.scaledacrossall`,
delta,
weight, survey_unit, season
)]
#add r-squared values
# Iterate through the best models and extract R-squared values
r_squared_values <-list()
for (i in 1:nrow(dd.dt.2.formatted)) {
summary_data <- summary(get.models(dd, subset = i)[[1]])
r_squared <- signif(summary_data$r.squared,2)
r_squared_values <- unlist(c(r_squared_values,r_squared))
}
dd.dt.2.formatted[,"R squared" := r_squared_values]
#empty data table
model_coef_se_fill <- data.table(Rank = as.numeric(), coef_name = as.character(),coef = as.numeric(), se = as.numeric())
for (i in 1:nrow(dd.dt.2.formatted)){
model_coef_se_single <- data.table(unlist(coefTable(dd,full = T)[[i]]))
model_coef_se_single[,coef_names := rownames(unlist(coefTable(dd,full = T)[[i]]))]
model_coef_se_single[,Rank := i]
colnames(model_coef_se_single) <- c("coef","se","df","coef_name","Rank")
#reduce to columns we need
model_coef_se_single <- model_coef_se_single[,.(Rank, coef_name, coef, se)]
model_coef_se_fill <- rbind(model_coef_se_fill,model_coef_se_single)
}
#format to merge with model rankings and averaged model
model_coef_se_fill[,coef_se := paste0(round(coef,3)," ± ",round(se,3))]
#delete extra columns
model_coef_se_fill <- model_coef_se_fill[,.(Rank,coef_name,coef_se)]
#long to wide
model_coef_se_fill.w <- dcast(model_coef_se_fill, formula = Rank ~ coef_name, value.var = c("coef_se"))
#merge model_coef_se_fill with dd.dt.2.formatted
model_coef_se_AIC <- dd.dt.2.formatted[model_coef_se_fill.w, on = "Rank"]
#model average all models with delta < 4 (8 models)
model_avg_delta4 <-model.avg(dd, subset = delta < 4, fit = T) #NB: The āsubsetā (or āconditionalā) average only averages over the models where the parameter appears. An alternative, the āfullā average assumes that a variable is included in every model, but in some models the corresponding coefficient (and its respective variance) is set to zero. Unlike the āsubset averageā, it does not have a tendency of biasing the value away from zero. The āfullā average is a type of shrinkage estimator, and for variables with a weak relationship to the response it is smaller than āsubsetā estimators., fit = T fits the component models again
model_avg_values <- as.data.table(coefTable(model_avg_delta4,fill = T)) # with SE
coef_names <- names(coef(model.avg(dd, subset = delta < 4)))
model_avg_values[,coef_name:=coef_names][,coef:=Estimate][,Estimate:=NULL][,df:=NULL][,se:= `Std. Error`][,`Std. Error` := NULL]
#new column with coef and SE
model_avg_values[,coef_se := paste0(round(coef,3)," ± ",round(se,3))]
#long to wide for model avg
model_avg.wide <- dcast(model_avg_values, formula = . ~ coef_name, value.var = c("coef_se"))
#add rank of "model avg" to table
model_avg.wide[,Rank := "Model avg"]
best_model_sbt_jaccard_table_formatted <- rbind(model_coef_se_AIC, model_avg.wide, fill = T)
#get rid of interaction coefficients
best_model_sbt_jaccard_table_formatted.r <- best_model_sbt_jaccard_table_formatted[,.(Rank,`(Intercept)`,
area_km.scaledacrossall,
depth_annual_avg.scaledacrossall,
depth_annual_range.scaledacrossall,
haul_id_count_annual.scaledacrossall,
latitude_annual_avg.scaledacrossall,
latitude_annual_range.scaledacrossall,
spp_count_annual.scaledacrossall,
summed_tonnes_scaled_byreg,
season,
survey_unit,
yearly_min_bypoint_avg.scaledacrossall,
`summed_tonnes_scaled_byreg:survey_unit`,
`survey_unit:yearly_min_bypoint_avg.scaledacrossall`,
`R squared`,
delta,
weight
)]
#round to 2 significant figures
#names of numeric columns
numeric_cols <- names(best_model_sbt_jaccard_table_formatted.r)[sapply(best_model_sbt_jaccard_table_formatted.r, is.numeric)]
# Apply signif() only to numeric columns
best_model_sbt_jaccard_table_formatted.r[, (numeric_cols) := lapply(.SD, function(x) if (is.numeric(x)) signif(x, digits = 2) else x), .SDcols = numeric_cols]
#change column names, in caption note that all variables are centered and scaled
colnames(best_model_sbt_jaccard_table_formatted.r) <- c(
"Rank",
"Intercept",
"Area", #scaled across all
"Average depth",
"Depth range",
"Number of tows",
"Average latitude",
"Latitude range",
"Species count",
"Relative catch", #scaled within region
"Season",
"Survey",
"Average minimum temperature",
"Survey * relative catch",
"Survey * avg min temperature",
"R squared",
paste0("\u0394"," AICc"),
paste0("\u03C9"))
#save as csv
fwrite(best_model_sbt_jaccard_table_formatted.r,here::here("output","best_model_sbt_jaccard_table_formatted.csv"))
###Predict dissimilarity across years using averaged model (model_avg_delta4)
dissimilarity_covariates_dredge.dt_predictions <- copy(dissimilarity_covariates_dredge.dt)
#allowing temp and fishing to vary (aka no changes)
dissimilarity_covariates_dredge.dt_predictions[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]] #full allows us to switch back to mixed effect models
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#What's the R^2 value?
r.squaredGLMM(lm(data = dissimilarity_covariates_dredge.dt_predictions, annual_dissimilarity_value ~ pred_dissim))
R2m R2c
[1,] 0.9445771 0.9445771
#R2 = 0.9507659
#constant temp in regions (aka take mean of temperature values within survey units so they are the same value within each year for a survey)
dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]
dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]]
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#What's the R^2 value?
r.squaredGLMM(lm(data = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg, annual_dissimilarity_value ~ pred_dissim))
R2m R2c
[1,] 0.941983 0.941983
#R2 = 0.94 (only fishing!!)
#constant fishing in regions (aka take mean of fishing values within survey units so they are the same value within each year for a survey)
dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)]
dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]]
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#What's the R^2 value?
r.squaredGLMM(lm(data = dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg, annual_dissimilarity_value ~ pred_dissim))
R2m R2c
[1,] 0.9336986 0.9336986
#R2 = 0.94 (only temperature!!)
#and then with consistent temp and fishing in regions (aka take mean of both min temp and fishing)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg <- copy(dissimilarity_covariates_dredge.dt)
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,summed_tonnes_scaled_byreg:=mean(summed_tonnes_scaled_byreg),.(survey_unit)][,yearly_min_bypoint_avg.scaledacrossall:=mean(yearly_min_bypoint_avg.scaledacrossall),.(survey_unit)]
dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,pred_dissim := predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg)[[1]]][,pred_se := predict(model_avg_delta4, se.fit = T, full = F)[[2]]]
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#What's the R^2 value?
r.squaredGLMM(lm(data = dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg, annual_dissimilarity_value ~ pred_dissim))
R2m R2c
[1,] 0.9314474 0.9314474
#R2 = 0.93 (only other variables!!! not temp or fishing)
#allowing temp and fishing to vary within regs (normal)
dissimilarity_covariates_dredge.dt_predictions[,pred_dissim :=
predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions)[[1]]][,pred_se :=
predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions)[[2]]]
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#allowing only fishing to vary within regions (with mean temp)
dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,pred_dissim :=
predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[1]]][,pred_se :=
predict(model_avg_delta4, se.fit = T, full = F, newdata = dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg)[[2]]]
Warning: argument 'full' ignoredWarning: argument 'full' ignored
#FOR COLOR TO MATCH
#sort color link by survey name season
#alphabetical order
color_link_alpha <- setorder(color_link, survey_unit)
#exclude surveys we don't include
color_link_alpha <- color_link_alpha[survey_unit %in% unique(dissimilarity_covariates_dredge.dt_predictions$survey_unit),]
color_alpha_order <- color_link_alpha[,hex]
label_alpha_order <- color_link_alpha[,Survey_Name_Season]
#maintain temp and fishing
#plot
predicted_values_temp_fishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions) +
geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
ylim(0,1.5) +
theme_classic() +
ggtitle("Average model predictions")
#average temp and fishing for each region
predicted_values_temp_fishing_meantempfishinginsurvey <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg) +
geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.3) +
scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
ylim(0,1.5) +
theme_classic() +
ggtitle("Average model predictions with mean\nsurvey temperature and fishing pressure")
#average temp and fishing across all regions
predicted_values_temp_fishing_meantempfishing <- ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing) +
geom_point(aes(x = year, y = annual_dissimilarity_value, color = survey_unit)) +
geom_line(aes(x = year, y = pred_dissim, color = survey_unit)) +
geom_ribbon(aes(x = year, ymin = pred_dissim-pred_se, ymax = pred_dissim+pred_se, fill = survey_unit), alpha = 0.1) +
scale_color_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
scale_fill_manual(values = color_alpha_order, labels = label_alpha_order, name = "Survey") +
labs(x = "Year",y = "Average annual total\nBray Curtis dissimilarity") +
ylim(0,1.5) +
theme_classic() +
ggtitle("Average model predictions with mean\novereall temperature and fishing pressure")
Error in ggplot(dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing) :
object 'dissimilarity_covariates_dredge.dt_predictions_consistenttempfishing' not found
Take dissimilarity values from random normal distribution for each year for each region, and then calculate slope (1000 times). Do this for: - Fishing and temperature vary interannually within surveys - Temperature is held constant (as mean over time series for a survey), but fishing varies (allows us to look at relative variance explained) - Fishing is held constant (as mean over time series for a survey), but temperature varies (allows us to look at relative variance explained) - Both fishing and temperature held constant (allows us to see role of other components of the model)
################################################################################
#full predictions
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table <- dissimilarity_covariates_dredge.dt_predictions[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
predicted_dissim_trends_rnormruns <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
table[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL TO EXTRACT SURVEY SPECIFIC TREND VALUES
jaccard_total_predicted_lm_singlerun <- lm(rnorm_pred ~ year*survey_unit,data = table)
model_coefs_reduced_predictions_singlerun <- data.table(summary(jaccard_total_predicted_lm_singlerun)$coefficients)
model_coefs_reduced_predictions_singlerun[,var := rownames(summary(jaccard_total_predicted_lm_singlerun)$coefficients)]
#limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
model_coefs_reduced_predictions_singlerun <- model_coefs_reduced_predictions_singlerun[c(2,34:64),]
#adjust survey unit name by deleting beginning of string
model_coefs_reduced_predictions_singlerun[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
#calculate interaction coefficients
AI_estimate <- model_coefs_reduced_predictions_singlerun[1,Estimate]
model_coefs_reduced_predictions_singlerun[1,estimate := AI_estimate]
model_coefs_reduced_predictions_singlerun[2:32,estimate := (AI_estimate + Estimate)]
predicted_dissim_trends_rnormruns <- rbind(predicted_dissim_trends_rnormruns, model_coefs_reduced_predictions_singlerun[,.(survey_unit, estimate)])
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212
[1] 213
[1] 214
[1] 215
[1] 216
[1] 217
[1] 218
[1] 219
[1] 220
[1] 221
[1] 222
[1] 223
[1] 224
[1] 225
[1] 226
[1] 227
[1] 228
[1] 229
[1] 230
[1] 231
[1] 232
[1] 233
[1] 234
[1] 235
[1] 236
[1] 237
[1] 238
[1] 239
[1] 240
[1] 241
[1] 242
[1] 243
[1] 244
[1] 245
[1] 246
[1] 247
[1] 248
[1] 249
[1] 250
[1] 251
[1] 252
[1] 253
[1] 254
[1] 255
[1] 256
[1] 257
[1] 258
[1] 259
[1] 260
[1] 261
[1] 262
[1] 263
[1] 264
[1] 265
[1] 266
[1] 267
[1] 268
[1] 269
[1] 270
[1] 271
[1] 272
[1] 273
[1] 274
[1] 275
[1] 276
[1] 277
[1] 278
[1] 279
[1] 280
[1] 281
[1] 282
[1] 283
[1] 284
[1] 285
[1] 286
[1] 287
[1] 288
[1] 289
[1] 290
[1] 291
[1] 292
[1] 293
[1] 294
[1] 295
[1] 296
[1] 297
[1] 298
[1] 299
[1] 300
[1] 301
[1] 302
[1] 303
[1] 304
[1] 305
[1] 306
[1] 307
[1] 308
[1] 309
[1] 310
[1] 311
[1] 312
[1] 313
[1] 314
[1] 315
[1] 316
[1] 317
[1] 318
[1] 319
[1] 320
[1] 321
[1] 322
[1] 323
[1] 324
[1] 325
[1] 326
[1] 327
[1] 328
[1] 329
[1] 330
[1] 331
[1] 332
[1] 333
[1] 334
[1] 335
[1] 336
[1] 337
[1] 338
[1] 339
[1] 340
[1] 341
[1] 342
[1] 343
[1] 344
[1] 345
[1] 346
[1] 347
[1] 348
[1] 349
[1] 350
[1] 351
[1] 352
[1] 353
[1] 354
[1] 355
[1] 356
[1] 357
[1] 358
[1] 359
[1] 360
[1] 361
[1] 362
[1] 363
[1] 364
[1] 365
[1] 366
[1] 367
[1] 368
[1] 369
[1] 370
[1] 371
[1] 372
[1] 373
[1] 374
[1] 375
[1] 376
[1] 377
[1] 378
[1] 379
[1] 380
[1] 381
[1] 382
[1] 383
[1] 384
[1] 385
[1] 386
[1] 387
[1] 388
[1] 389
[1] 390
[1] 391
[1] 392
[1] 393
[1] 394
[1] 395
[1] 396
[1] 397
[1] 398
[1] 399
[1] 400
[1] 401
[1] 402
[1] 403
[1] 404
[1] 405
[1] 406
[1] 407
[1] 408
[1] 409
[1] 410
[1] 411
[1] 412
[1] 413
[1] 414
[1] 415
[1] 416
[1] 417
[1] 418
[1] 419
[1] 420
[1] 421
[1] 422
[1] 423
[1] 424
[1] 425
[1] 426
[1] 427
[1] 428
[1] 429
[1] 430
[1] 431
[1] 432
[1] 433
[1] 434
[1] 435
[1] 436
[1] 437
[1] 438
[1] 439
[1] 440
[1] 441
[1] 442
[1] 443
[1] 444
[1] 445
[1] 446
[1] 447
[1] 448
[1] 449
[1] 450
[1] 451
[1] 452
[1] 453
[1] 454
[1] 455
[1] 456
[1] 457
[1] 458
[1] 459
[1] 460
[1] 461
[1] 462
[1] 463
[1] 464
[1] 465
[1] 466
[1] 467
[1] 468
[1] 469
[1] 470
[1] 471
[1] 472
[1] 473
[1] 474
[1] 475
[1] 476
[1] 477
[1] 478
[1] 479
[1] 480
[1] 481
[1] 482
[1] 483
[1] 484
[1] 485
[1] 486
[1] 487
[1] 488
[1] 489
[1] 490
[1] 491
[1] 492
[1] 493
[1] 494
[1] 495
[1] 496
[1] 497
[1] 498
[1] 499
[1] 500
[1] 501
[1] 502
[1] 503
[1] 504
[1] 505
[1] 506
[1] 507
[1] 508
[1] 509
[1] 510
[1] 511
[1] 512
[1] 513
[1] 514
[1] 515
[1] 516
[1] 517
[1] 518
[1] 519
[1] 520
[1] 521
[1] 522
[1] 523
[1] 524
[1] 525
[1] 526
[1] 527
[1] 528
[1] 529
[1] 530
[1] 531
[1] 532
[1] 533
[1] 534
[1] 535
[1] 536
[1] 537
[1] 538
[1] 539
[1] 540
[1] 541
[1] 542
[1] 543
[1] 544
[1] 545
[1] 546
[1] 547
[1] 548
[1] 549
[1] 550
[1] 551
[1] 552
[1] 553
[1] 554
[1] 555
[1] 556
[1] 557
[1] 558
[1] 559
[1] 560
[1] 561
[1] 562
[1] 563
[1] 564
[1] 565
[1] 566
[1] 567
[1] 568
[1] 569
[1] 570
[1] 571
[1] 572
[1] 573
[1] 574
[1] 575
[1] 576
[1] 577
[1] 578
[1] 579
[1] 580
[1] 581
[1] 582
[1] 583
[1] 584
[1] 585
[1] 586
[1] 587
[1] 588
[1] 589
[1] 590
[1] 591
[1] 592
[1] 593
[1] 594
[1] 595
[1] 596
[1] 597
[1] 598
[1] 599
[1] 600
[1] 601
[1] 602
[1] 603
[1] 604
[1] 605
[1] 606
[1] 607
[1] 608
[1] 609
[1] 610
[1] 611
[1] 612
[1] 613
[1] 614
[1] 615
[1] 616
[1] 617
[1] 618
[1] 619
[1] 620
[1] 621
[1] 622
[1] 623
[1] 624
[1] 625
[1] 626
[1] 627
[1] 628
[1] 629
[1] 630
[1] 631
[1] 632
[1] 633
[1] 634
[1] 635
[1] 636
[1] 637
[1] 638
[1] 639
[1] 640
[1] 641
[1] 642
[1] 643
[1] 644
[1] 645
[1] 646
[1] 647
[1] 648
[1] 649
[1] 650
[1] 651
[1] 652
[1] 653
[1] 654
[1] 655
[1] 656
[1] 657
[1] 658
[1] 659
[1] 660
[1] 661
[1] 662
[1] 663
[1] 664
[1] 665
[1] 666
[1] 667
[1] 668
[1] 669
[1] 670
[1] 671
[1] 672
[1] 673
[1] 674
[1] 675
[1] 676
[1] 677
[1] 678
[1] 679
[1] 680
[1] 681
[1] 682
[1] 683
[1] 684
[1] 685
[1] 686
[1] 687
[1] 688
[1] 689
[1] 690
[1] 691
[1] 692
[1] 693
[1] 694
[1] 695
[1] 696
[1] 697
[1] 698
[1] 699
[1] 700
[1] 701
[1] 702
[1] 703
[1] 704
[1] 705
[1] 706
[1] 707
[1] 708
[1] 709
[1] 710
[1] 711
[1] 712
[1] 713
[1] 714
[1] 715
[1] 716
[1] 717
[1] 718
[1] 719
[1] 720
[1] 721
[1] 722
[1] 723
[1] 724
[1] 725
[1] 726
[1] 727
[1] 728
[1] 729
[1] 730
[1] 731
[1] 732
[1] 733
[1] 734
[1] 735
[1] 736
[1] 737
[1] 738
[1] 739
[1] 740
[1] 741
[1] 742
[1] 743
[1] 744
[1] 745
[1] 746
[1] 747
[1] 748
[1] 749
[1] 750
[1] 751
[1] 752
[1] 753
[1] 754
[1] 755
[1] 756
[1] 757
[1] 758
[1] 759
[1] 760
[1] 761
[1] 762
[1] 763
[1] 764
[1] 765
[1] 766
[1] 767
[1] 768
[1] 769
[1] 770
[1] 771
[1] 772
[1] 773
[1] 774
[1] 775
[1] 776
[1] 777
[1] 778
[1] 779
[1] 780
[1] 781
[1] 782
[1] 783
[1] 784
[1] 785
[1] 786
[1] 787
[1] 788
[1] 789
[1] 790
[1] 791
[1] 792
[1] 793
[1] 794
[1] 795
[1] 796
[1] 797
[1] 798
[1] 799
[1] 800
[1] 801
[1] 802
[1] 803
[1] 804
[1] 805
[1] 806
[1] 807
[1] 808
[1] 809
[1] 810
[1] 811
[1] 812
[1] 813
[1] 814
[1] 815
[1] 816
[1] 817
[1] 818
[1] 819
[1] 820
[1] 821
[1] 822
[1] 823
[1] 824
[1] 825
[1] 826
[1] 827
[1] 828
[1] 829
[1] 830
[1] 831
[1] 832
[1] 833
[1] 834
[1] 835
[1] 836
[1] 837
[1] 838
[1] 839
[1] 840
[1] 841
[1] 842
[1] 843
[1] 844
[1] 845
[1] 846
[1] 847
[1] 848
[1] 849
[1] 850
[1] 851
[1] 852
[1] 853
[1] 854
[1] 855
[1] 856
[1] 857
[1] 858
[1] 859
[1] 860
[1] 861
[1] 862
[1] 863
[1] 864
[1] 865
[1] 866
[1] 867
[1] 868
[1] 869
[1] 870
[1] 871
[1] 872
[1] 873
[1] 874
[1] 875
[1] 876
[1] 877
[1] 878
[1] 879
[1] 880
[1] 881
[1] 882
[1] 883
[1] 884
[1] 885
[1] 886
[1] 887
[1] 888
[1] 889
[1] 890
[1] 891
[1] 892
[1] 893
[1] 894
[1] 895
[1] 896
[1] 897
[1] 898
[1] 899
[1] 900
[1] 901
[1] 902
[1] 903
[1] 904
[1] 905
[1] 906
[1] 907
[1] 908
[1] 909
[1] 910
[1] 911
[1] 912
[1] 913
[1] 914
[1] 915
[1] 916
[1] 917
[1] 918
[1] 919
[1] 920
[1] 921
[1] 922
[1] 923
[1] 924
[1] 925
[1] 926
[1] 927
[1] 928
[1] 929
[1] 930
[1] 931
[1] 932
[1] 933
[1] 934
[1] 935
[1] 936
[1] 937
[1] 938
[1] 939
[1] 940
[1] 941
[1] 942
[1] 943
[1] 944
[1] 945
[1] 946
[1] 947
[1] 948
[1] 949
[1] 950
[1] 951
[1] 952
[1] 953
[1] 954
[1] 955
[1] 956
[1] 957
[1] 958
[1] 959
[1] 960
[1] 961
[1] 962
[1] 963
[1] 964
[1] 965
[1] 966
[1] 967
[1] 968
[1] 969
[1] 970
[1] 971
[1] 972
[1] 973
[1] 974
[1] 975
[1] 976
[1] 977
[1] 978
[1] 979
[1] 980
[1] 981
[1] 982
[1] 983
[1] 984
[1] 985
[1] 986
[1] 987
[1] 988
[1] 989
[1] 990
[1] 991
[1] 992
[1] 993
[1] 994
[1] 995
[1] 996
[1] 997
[1] 998
[1] 999
[1] 1000
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]
predicted_dissim_trends_rnormruns.summary <- unique(predicted_dissim_trends_rnormruns[,.(survey_unit, mean_dissim_coef, sd_dissim)])
predicted_dissim_trends_rnormruns.summary[,pred_type := "full"]
################################################################################
#predictions with temp held constant and fishing still varying from year to year
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttemp <- dissimilarity_covariates_dredge.dt_predictions_consistenttempinreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
predicted_dissim_trends_rnormruns_constanttemp <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
table_constanttemp[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
jaccard_total_predicted_lm_singlerun_constanttemp <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttemp)
model_coefs_reduced_predictions_singlerun_constanttemp <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)
model_coefs_reduced_predictions_singlerun_constanttemp[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttemp)$coefficients)]
#limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
model_coefs_reduced_predictions_singlerun_constanttemp <- model_coefs_reduced_predictions_singlerun_constanttemp[c(2,34:64),]
#adjust survey unit name by deleting beginning of string
model_coefs_reduced_predictions_singlerun_constanttemp[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
#calculate interaction coefficients
AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttemp[1,Estimate]
model_coefs_reduced_predictions_singlerun_constanttemp[1,estimate := AI_estimate]
model_coefs_reduced_predictions_singlerun_constanttemp[2:32,estimate := (AI_estimate + Estimate)]
predicted_dissim_trends_rnormruns_constanttemp <- rbind(predicted_dissim_trends_rnormruns_constanttemp, model_coefs_reduced_predictions_singlerun_constanttemp[,.(survey_unit, estimate)])
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212
[1] 213
[1] 214
[1] 215
[1] 216
[1] 217
[1] 218
[1] 219
[1] 220
[1] 221
[1] 222
[1] 223
[1] 224
[1] 225
[1] 226
[1] 227
[1] 228
[1] 229
[1] 230
[1] 231
[1] 232
[1] 233
[1] 234
[1] 235
[1] 236
[1] 237
[1] 238
[1] 239
[1] 240
[1] 241
[1] 242
[1] 243
[1] 244
[1] 245
[1] 246
[1] 247
[1] 248
[1] 249
[1] 250
[1] 251
[1] 252
[1] 253
[1] 254
[1] 255
[1] 256
[1] 257
[1] 258
[1] 259
[1] 260
[1] 261
[1] 262
[1] 263
[1] 264
[1] 265
[1] 266
[1] 267
[1] 268
[1] 269
[1] 270
[1] 271
[1] 272
[1] 273
[1] 274
[1] 275
[1] 276
[1] 277
[1] 278
[1] 279
[1] 280
[1] 281
[1] 282
[1] 283
[1] 284
[1] 285
[1] 286
[1] 287
[1] 288
[1] 289
[1] 290
[1] 291
[1] 292
[1] 293
[1] 294
[1] 295
[1] 296
[1] 297
[1] 298
[1] 299
[1] 300
[1] 301
[1] 302
[1] 303
[1] 304
[1] 305
[1] 306
[1] 307
[1] 308
[1] 309
[1] 310
[1] 311
[1] 312
[1] 313
[1] 314
[1] 315
[1] 316
[1] 317
[1] 318
[1] 319
[1] 320
[1] 321
[1] 322
[1] 323
[1] 324
[1] 325
[1] 326
[1] 327
[1] 328
[1] 329
[1] 330
[1] 331
[1] 332
[1] 333
[1] 334
[1] 335
[1] 336
[1] 337
[1] 338
[1] 339
[1] 340
[1] 341
[1] 342
[1] 343
[1] 344
[1] 345
[1] 346
[1] 347
[1] 348
[1] 349
[1] 350
[1] 351
[1] 352
[1] 353
[1] 354
[1] 355
[1] 356
[1] 357
[1] 358
[1] 359
[1] 360
[1] 361
[1] 362
[1] 363
[1] 364
[1] 365
[1] 366
[1] 367
[1] 368
[1] 369
[1] 370
[1] 371
[1] 372
[1] 373
[1] 374
[1] 375
[1] 376
[1] 377
[1] 378
[1] 379
[1] 380
[1] 381
[1] 382
[1] 383
[1] 384
[1] 385
[1] 386
[1] 387
[1] 388
[1] 389
[1] 390
[1] 391
[1] 392
[1] 393
[1] 394
[1] 395
[1] 396
[1] 397
[1] 398
[1] 399
[1] 400
[1] 401
[1] 402
[1] 403
[1] 404
[1] 405
[1] 406
[1] 407
[1] 408
[1] 409
[1] 410
[1] 411
[1] 412
[1] 413
[1] 414
[1] 415
[1] 416
[1] 417
[1] 418
[1] 419
[1] 420
[1] 421
[1] 422
[1] 423
[1] 424
[1] 425
[1] 426
[1] 427
[1] 428
[1] 429
[1] 430
[1] 431
[1] 432
[1] 433
[1] 434
[1] 435
[1] 436
[1] 437
[1] 438
[1] 439
[1] 440
[1] 441
[1] 442
[1] 443
[1] 444
[1] 445
[1] 446
[1] 447
[1] 448
[1] 449
[1] 450
[1] 451
[1] 452
[1] 453
[1] 454
[1] 455
[1] 456
[1] 457
[1] 458
[1] 459
[1] 460
[1] 461
[1] 462
[1] 463
[1] 464
[1] 465
[1] 466
[1] 467
[1] 468
[1] 469
[1] 470
[1] 471
[1] 472
[1] 473
[1] 474
[1] 475
[1] 476
[1] 477
[1] 478
[1] 479
[1] 480
[1] 481
[1] 482
[1] 483
[1] 484
[1] 485
[1] 486
[1] 487
[1] 488
[1] 489
[1] 490
[1] 491
[1] 492
[1] 493
[1] 494
[1] 495
[1] 496
[1] 497
[1] 498
[1] 499
[1] 500
[1] 501
[1] 502
[1] 503
[1] 504
[1] 505
[1] 506
[1] 507
[1] 508
[1] 509
[1] 510
[1] 511
[1] 512
[1] 513
[1] 514
[1] 515
[1] 516
[1] 517
[1] 518
[1] 519
[1] 520
[1] 521
[1] 522
[1] 523
[1] 524
[1] 525
[1] 526
[1] 527
[1] 528
[1] 529
[1] 530
[1] 531
[1] 532
[1] 533
[1] 534
[1] 535
[1] 536
[1] 537
[1] 538
[1] 539
[1] 540
[1] 541
[1] 542
[1] 543
[1] 544
[1] 545
[1] 546
[1] 547
[1] 548
[1] 549
[1] 550
[1] 551
[1] 552
[1] 553
[1] 554
[1] 555
[1] 556
[1] 557
[1] 558
[1] 559
[1] 560
[1] 561
[1] 562
[1] 563
[1] 564
[1] 565
[1] 566
[1] 567
[1] 568
[1] 569
[1] 570
[1] 571
[1] 572
[1] 573
[1] 574
[1] 575
[1] 576
[1] 577
[1] 578
[1] 579
[1] 580
[1] 581
[1] 582
[1] 583
[1] 584
[1] 585
[1] 586
[1] 587
[1] 588
[1] 589
[1] 590
[1] 591
[1] 592
[1] 593
[1] 594
[1] 595
[1] 596
[1] 597
[1] 598
[1] 599
[1] 600
[1] 601
[1] 602
[1] 603
[1] 604
[1] 605
[1] 606
[1] 607
[1] 608
[1] 609
[1] 610
[1] 611
[1] 612
[1] 613
[1] 614
[1] 615
[1] 616
[1] 617
[1] 618
[1] 619
[1] 620
[1] 621
[1] 622
[1] 623
[1] 624
[1] 625
[1] 626
[1] 627
[1] 628
[1] 629
[1] 630
[1] 631
[1] 632
[1] 633
[1] 634
[1] 635
[1] 636
[1] 637
[1] 638
[1] 639
[1] 640
[1] 641
[1] 642
[1] 643
[1] 644
[1] 645
[1] 646
[1] 647
[1] 648
[1] 649
[1] 650
[1] 651
[1] 652
[1] 653
[1] 654
[1] 655
[1] 656
[1] 657
[1] 658
[1] 659
[1] 660
[1] 661
[1] 662
[1] 663
[1] 664
[1] 665
[1] 666
[1] 667
[1] 668
[1] 669
[1] 670
[1] 671
[1] 672
[1] 673
[1] 674
[1] 675
[1] 676
[1] 677
[1] 678
[1] 679
[1] 680
[1] 681
[1] 682
[1] 683
[1] 684
[1] 685
[1] 686
[1] 687
[1] 688
[1] 689
[1] 690
[1] 691
[1] 692
[1] 693
[1] 694
[1] 695
[1] 696
[1] 697
[1] 698
[1] 699
[1] 700
[1] 701
[1] 702
[1] 703
[1] 704
[1] 705
[1] 706
[1] 707
[1] 708
[1] 709
[1] 710
[1] 711
[1] 712
[1] 713
[1] 714
[1] 715
[1] 716
[1] 717
[1] 718
[1] 719
[1] 720
[1] 721
[1] 722
[1] 723
[1] 724
[1] 725
[1] 726
[1] 727
[1] 728
[1] 729
[1] 730
[1] 731
[1] 732
[1] 733
[1] 734
[1] 735
[1] 736
[1] 737
[1] 738
[1] 739
[1] 740
[1] 741
[1] 742
[1] 743
[1] 744
[1] 745
[1] 746
[1] 747
[1] 748
[1] 749
[1] 750
[1] 751
[1] 752
[1] 753
[1] 754
[1] 755
[1] 756
[1] 757
[1] 758
[1] 759
[1] 760
[1] 761
[1] 762
[1] 763
[1] 764
[1] 765
[1] 766
[1] 767
[1] 768
[1] 769
[1] 770
[1] 771
[1] 772
[1] 773
[1] 774
[1] 775
[1] 776
[1] 777
[1] 778
[1] 779
[1] 780
[1] 781
[1] 782
[1] 783
[1] 784
[1] 785
[1] 786
[1] 787
[1] 788
[1] 789
[1] 790
[1] 791
[1] 792
[1] 793
[1] 794
[1] 795
[1] 796
[1] 797
[1] 798
[1] 799
[1] 800
[1] 801
[1] 802
[1] 803
[1] 804
[1] 805
[1] 806
[1] 807
[1] 808
[1] 809
[1] 810
[1] 811
[1] 812
[1] 813
[1] 814
[1] 815
[1] 816
[1] 817
[1] 818
[1] 819
[1] 820
[1] 821
[1] 822
[1] 823
[1] 824
[1] 825
[1] 826
[1] 827
[1] 828
[1] 829
[1] 830
[1] 831
[1] 832
[1] 833
[1] 834
[1] 835
[1] 836
[1] 837
[1] 838
[1] 839
[1] 840
[1] 841
[1] 842
[1] 843
[1] 844
[1] 845
[1] 846
[1] 847
[1] 848
[1] 849
[1] 850
[1] 851
[1] 852
[1] 853
[1] 854
[1] 855
[1] 856
[1] 857
[1] 858
[1] 859
[1] 860
[1] 861
[1] 862
[1] 863
[1] 864
[1] 865
[1] 866
[1] 867
[1] 868
[1] 869
[1] 870
[1] 871
[1] 872
[1] 873
[1] 874
[1] 875
[1] 876
[1] 877
[1] 878
[1] 879
[1] 880
[1] 881
[1] 882
[1] 883
[1] 884
[1] 885
[1] 886
[1] 887
[1] 888
[1] 889
[1] 890
[1] 891
[1] 892
[1] 893
[1] 894
[1] 895
[1] 896
[1] 897
[1] 898
[1] 899
[1] 900
[1] 901
[1] 902
[1] 903
[1] 904
[1] 905
[1] 906
[1] 907
[1] 908
[1] 909
[1] 910
[1] 911
[1] 912
[1] 913
[1] 914
[1] 915
[1] 916
[1] 917
[1] 918
[1] 919
[1] 920
[1] 921
[1] 922
[1] 923
[1] 924
[1] 925
[1] 926
[1] 927
[1] 928
[1] 929
[1] 930
[1] 931
[1] 932
[1] 933
[1] 934
[1] 935
[1] 936
[1] 937
[1] 938
[1] 939
[1] 940
[1] 941
[1] 942
[1] 943
[1] 944
[1] 945
[1] 946
[1] 947
[1] 948
[1] 949
[1] 950
[1] 951
[1] 952
[1] 953
[1] 954
[1] 955
[1] 956
[1] 957
[1] 958
[1] 959
[1] 960
[1] 961
[1] 962
[1] 963
[1] 964
[1] 965
[1] 966
[1] 967
[1] 968
[1] 969
[1] 970
[1] 971
[1] 972
[1] 973
[1] 974
[1] 975
[1] 976
[1] 977
[1] 978
[1] 979
[1] 980
[1] 981
[1] 982
[1] 983
[1] 984
[1] 985
[1] 986
[1] 987
[1] 988
[1] 989
[1] 990
[1] 991
[1] 992
[1] 993
[1] 994
[1] 995
[1] 996
[1] 997
[1] 998
[1] 999
[1] 1000
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttemp[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]
predicted_dissim_trends_rnormrunsconstant_temp.summary <- unique(predicted_dissim_trends_rnormruns_constanttemp[,.(survey_unit, mean_dissim_coef, sd_dissim)])
predicted_dissim_trends_rnormrunsconstant_temp.summary[,pred_type := "temp_constant"]
predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_temp.summary)
################################################################################
#predictions with fishing held constant (and temperature varying)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constantfishing <- dissimilarity_covariates_dredge.dt_predictions_consistentfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
predicted_dissim_trends_rnormruns_constantfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
table_constantfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
jaccard_total_predicted_lm_singlerun_constantfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constantfishing)
model_coefs_reduced_predictions_singlerun_constantfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)
model_coefs_reduced_predictions_singlerun_constantfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constantfishing)$coefficients)]
#limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
model_coefs_reduced_predictions_singlerun_constantfishing <- model_coefs_reduced_predictions_singlerun_constantfishing[c(2,34:64),]
#adjust survey unit name by deleting beginning of string
model_coefs_reduced_predictions_singlerun_constantfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
#calculate interaction coefficients
AI_estimate <- model_coefs_reduced_predictions_singlerun_constantfishing[1,Estimate]
model_coefs_reduced_predictions_singlerun_constantfishing[1,estimate := AI_estimate]
model_coefs_reduced_predictions_singlerun_constantfishing[2:32,estimate := (AI_estimate + Estimate)]
predicted_dissim_trends_rnormruns_constantfishing <- rbind(predicted_dissim_trends_rnormruns_constantfishing, model_coefs_reduced_predictions_singlerun_constantfishing[,.(survey_unit, estimate)])
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212
[1] 213
[1] 214
[1] 215
[1] 216
[1] 217
[1] 218
[1] 219
[1] 220
[1] 221
[1] 222
[1] 223
[1] 224
[1] 225
[1] 226
[1] 227
[1] 228
[1] 229
[1] 230
[1] 231
[1] 232
[1] 233
[1] 234
[1] 235
[1] 236
[1] 237
[1] 238
[1] 239
[1] 240
[1] 241
[1] 242
[1] 243
[1] 244
[1] 245
[1] 246
[1] 247
[1] 248
[1] 249
[1] 250
[1] 251
[1] 252
[1] 253
[1] 254
[1] 255
[1] 256
[1] 257
[1] 258
[1] 259
[1] 260
[1] 261
[1] 262
[1] 263
[1] 264
[1] 265
[1] 266
[1] 267
[1] 268
[1] 269
[1] 270
[1] 271
[1] 272
[1] 273
[1] 274
[1] 275
[1] 276
[1] 277
[1] 278
[1] 279
[1] 280
[1] 281
[1] 282
[1] 283
[1] 284
[1] 285
[1] 286
[1] 287
[1] 288
[1] 289
[1] 290
[1] 291
[1] 292
[1] 293
[1] 294
[1] 295
[1] 296
[1] 297
[1] 298
[1] 299
[1] 300
[1] 301
[1] 302
[1] 303
[1] 304
[1] 305
[1] 306
[1] 307
[1] 308
[1] 309
[1] 310
[1] 311
[1] 312
[1] 313
[1] 314
[1] 315
[1] 316
[1] 317
[1] 318
[1] 319
[1] 320
[1] 321
[1] 322
[1] 323
[1] 324
[1] 325
[1] 326
[1] 327
[1] 328
[1] 329
[1] 330
[1] 331
[1] 332
[1] 333
[1] 334
[1] 335
[1] 336
[1] 337
[1] 338
[1] 339
[1] 340
[1] 341
[1] 342
[1] 343
[1] 344
[1] 345
[1] 346
[1] 347
[1] 348
[1] 349
[1] 350
[1] 351
[1] 352
[1] 353
[1] 354
[1] 355
[1] 356
[1] 357
[1] 358
[1] 359
[1] 360
[1] 361
[1] 362
[1] 363
[1] 364
[1] 365
[1] 366
[1] 367
[1] 368
[1] 369
[1] 370
[1] 371
[1] 372
[1] 373
[1] 374
[1] 375
[1] 376
[1] 377
[1] 378
[1] 379
[1] 380
[1] 381
[1] 382
[1] 383
[1] 384
[1] 385
[1] 386
[1] 387
[1] 388
[1] 389
[1] 390
[1] 391
[1] 392
[1] 393
[1] 394
[1] 395
[1] 396
[1] 397
[1] 398
[1] 399
[1] 400
[1] 401
[1] 402
[1] 403
[1] 404
[1] 405
[1] 406
[1] 407
[1] 408
[1] 409
[1] 410
[1] 411
[1] 412
[1] 413
[1] 414
[1] 415
[1] 416
[1] 417
[1] 418
[1] 419
[1] 420
[1] 421
[1] 422
[1] 423
[1] 424
[1] 425
[1] 426
[1] 427
[1] 428
[1] 429
[1] 430
[1] 431
[1] 432
[1] 433
[1] 434
[1] 435
[1] 436
[1] 437
[1] 438
[1] 439
[1] 440
[1] 441
[1] 442
[1] 443
[1] 444
[1] 445
[1] 446
[1] 447
[1] 448
[1] 449
[1] 450
[1] 451
[1] 452
[1] 453
[1] 454
[1] 455
[1] 456
[1] 457
[1] 458
[1] 459
[1] 460
[1] 461
[1] 462
[1] 463
[1] 464
[1] 465
[1] 466
[1] 467
[1] 468
[1] 469
[1] 470
[1] 471
[1] 472
[1] 473
[1] 474
[1] 475
[1] 476
[1] 477
[1] 478
[1] 479
[1] 480
[1] 481
[1] 482
[1] 483
[1] 484
[1] 485
[1] 486
[1] 487
[1] 488
[1] 489
[1] 490
[1] 491
[1] 492
[1] 493
[1] 494
[1] 495
[1] 496
[1] 497
[1] 498
[1] 499
[1] 500
[1] 501
[1] 502
[1] 503
[1] 504
[1] 505
[1] 506
[1] 507
[1] 508
[1] 509
[1] 510
[1] 511
[1] 512
[1] 513
[1] 514
[1] 515
[1] 516
[1] 517
[1] 518
[1] 519
[1] 520
[1] 521
[1] 522
[1] 523
[1] 524
[1] 525
[1] 526
[1] 527
[1] 528
[1] 529
[1] 530
[1] 531
[1] 532
[1] 533
[1] 534
[1] 535
[1] 536
[1] 537
[1] 538
[1] 539
[1] 540
[1] 541
[1] 542
[1] 543
[1] 544
[1] 545
[1] 546
[1] 547
[1] 548
[1] 549
[1] 550
[1] 551
[1] 552
[1] 553
[1] 554
[1] 555
[1] 556
[1] 557
[1] 558
[1] 559
[1] 560
[1] 561
[1] 562
[1] 563
[1] 564
[1] 565
[1] 566
[1] 567
[1] 568
[1] 569
[1] 570
[1] 571
[1] 572
[1] 573
[1] 574
[1] 575
[1] 576
[1] 577
[1] 578
[1] 579
[1] 580
[1] 581
[1] 582
[1] 583
[1] 584
[1] 585
[1] 586
[1] 587
[1] 588
[1] 589
[1] 590
[1] 591
[1] 592
[1] 593
[1] 594
[1] 595
[1] 596
[1] 597
[1] 598
[1] 599
[1] 600
[1] 601
[1] 602
[1] 603
[1] 604
[1] 605
[1] 606
[1] 607
[1] 608
[1] 609
[1] 610
[1] 611
[1] 612
[1] 613
[1] 614
[1] 615
[1] 616
[1] 617
[1] 618
[1] 619
[1] 620
[1] 621
[1] 622
[1] 623
[1] 624
[1] 625
[1] 626
[1] 627
[1] 628
[1] 629
[1] 630
[1] 631
[1] 632
[1] 633
[1] 634
[1] 635
[1] 636
[1] 637
[1] 638
[1] 639
[1] 640
[1] 641
[1] 642
[1] 643
[1] 644
[1] 645
[1] 646
[1] 647
[1] 648
[1] 649
[1] 650
[1] 651
[1] 652
[1] 653
[1] 654
[1] 655
[1] 656
[1] 657
[1] 658
[1] 659
[1] 660
[1] 661
[1] 662
[1] 663
[1] 664
[1] 665
[1] 666
[1] 667
[1] 668
[1] 669
[1] 670
[1] 671
[1] 672
[1] 673
[1] 674
[1] 675
[1] 676
[1] 677
[1] 678
[1] 679
[1] 680
[1] 681
[1] 682
[1] 683
[1] 684
[1] 685
[1] 686
[1] 687
[1] 688
[1] 689
[1] 690
[1] 691
[1] 692
[1] 693
[1] 694
[1] 695
[1] 696
[1] 697
[1] 698
[1] 699
[1] 700
[1] 701
[1] 702
[1] 703
[1] 704
[1] 705
[1] 706
[1] 707
[1] 708
[1] 709
[1] 710
[1] 711
[1] 712
[1] 713
[1] 714
[1] 715
[1] 716
[1] 717
[1] 718
[1] 719
[1] 720
[1] 721
[1] 722
[1] 723
[1] 724
[1] 725
[1] 726
[1] 727
[1] 728
[1] 729
[1] 730
[1] 731
[1] 732
[1] 733
[1] 734
[1] 735
[1] 736
[1] 737
[1] 738
[1] 739
[1] 740
[1] 741
[1] 742
[1] 743
[1] 744
[1] 745
[1] 746
[1] 747
[1] 748
[1] 749
[1] 750
[1] 751
[1] 752
[1] 753
[1] 754
[1] 755
[1] 756
[1] 757
[1] 758
[1] 759
[1] 760
[1] 761
[1] 762
[1] 763
[1] 764
[1] 765
[1] 766
[1] 767
[1] 768
[1] 769
[1] 770
[1] 771
[1] 772
[1] 773
[1] 774
[1] 775
[1] 776
[1] 777
[1] 778
[1] 779
[1] 780
[1] 781
[1] 782
[1] 783
[1] 784
[1] 785
[1] 786
[1] 787
[1] 788
[1] 789
[1] 790
[1] 791
[1] 792
[1] 793
[1] 794
[1] 795
[1] 796
[1] 797
[1] 798
[1] 799
[1] 800
[1] 801
[1] 802
[1] 803
[1] 804
[1] 805
[1] 806
[1] 807
[1] 808
[1] 809
[1] 810
[1] 811
[1] 812
[1] 813
[1] 814
[1] 815
[1] 816
[1] 817
[1] 818
[1] 819
[1] 820
[1] 821
[1] 822
[1] 823
[1] 824
[1] 825
[1] 826
[1] 827
[1] 828
[1] 829
[1] 830
[1] 831
[1] 832
[1] 833
[1] 834
[1] 835
[1] 836
[1] 837
[1] 838
[1] 839
[1] 840
[1] 841
[1] 842
[1] 843
[1] 844
[1] 845
[1] 846
[1] 847
[1] 848
[1] 849
[1] 850
[1] 851
[1] 852
[1] 853
[1] 854
[1] 855
[1] 856
[1] 857
[1] 858
[1] 859
[1] 860
[1] 861
[1] 862
[1] 863
[1] 864
[1] 865
[1] 866
[1] 867
[1] 868
[1] 869
[1] 870
[1] 871
[1] 872
[1] 873
[1] 874
[1] 875
[1] 876
[1] 877
[1] 878
[1] 879
[1] 880
[1] 881
[1] 882
[1] 883
[1] 884
[1] 885
[1] 886
[1] 887
[1] 888
[1] 889
[1] 890
[1] 891
[1] 892
[1] 893
[1] 894
[1] 895
[1] 896
[1] 897
[1] 898
[1] 899
[1] 900
[1] 901
[1] 902
[1] 903
[1] 904
[1] 905
[1] 906
[1] 907
[1] 908
[1] 909
[1] 910
[1] 911
[1] 912
[1] 913
[1] 914
[1] 915
[1] 916
[1] 917
[1] 918
[1] 919
[1] 920
[1] 921
[1] 922
[1] 923
[1] 924
[1] 925
[1] 926
[1] 927
[1] 928
[1] 929
[1] 930
[1] 931
[1] 932
[1] 933
[1] 934
[1] 935
[1] 936
[1] 937
[1] 938
[1] 939
[1] 940
[1] 941
[1] 942
[1] 943
[1] 944
[1] 945
[1] 946
[1] 947
[1] 948
[1] 949
[1] 950
[1] 951
[1] 952
[1] 953
[1] 954
[1] 955
[1] 956
[1] 957
[1] 958
[1] 959
[1] 960
[1] 961
[1] 962
[1] 963
[1] 964
[1] 965
[1] 966
[1] 967
[1] 968
[1] 969
[1] 970
[1] 971
[1] 972
[1] 973
[1] 974
[1] 975
[1] 976
[1] 977
[1] 978
[1] 979
[1] 980
[1] 981
[1] 982
[1] 983
[1] 984
[1] 985
[1] 986
[1] 987
[1] 988
[1] 989
[1] 990
[1] 991
[1] 992
[1] 993
[1] 994
[1] 995
[1] 996
[1] 997
[1] 998
[1] 999
[1] 1000
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constantfishing[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]
predicted_dissim_trends_rnormrunsconstant_fishing.summary <- unique(predicted_dissim_trends_rnormruns_constantfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])
predicted_dissim_trends_rnormrunsconstant_fishing.summary[,pred_type := "fishing_constant"]
predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_fishing.summary)
################################################################################
#predictions with both fishing and temperature held constant (variability goes to other factors we don't account for)
################################################################################
#table with predicted dissimilarity values and standard error of all predicted dissimilarity values (by year)
table_constanttempfishing <- dissimilarity_covariates_dredge.dt_predictions_consistenttempfishinginreg[,.(survey_unit, pred_dissim, pred_se, year)]
#0) make datatable to populate
predicted_dissim_trends_rnormruns_constanttempfishing <- data.table()
#1) NEW PREDICTED VALUES FROM DISTRIBUTION
for (i in 1:1000){
table_constanttempfishing[,rnorm_pred := rnorm(1, mean = pred_dissim, sd = pred_se),.(year, survey_unit)]
#2) CALCULATE LINEAR MODEL FOR SLOPE VALUES
jaccard_total_predicted_lm_singlerun_constanttempfishing <- lm(rnorm_pred ~ year*survey_unit,data = table_constanttempfishing)
model_coefs_reduced_predictions_singlerun_constanttempfishing <- data.table(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)
model_coefs_reduced_predictions_singlerun_constanttempfishing[,var := rownames(summary(jaccard_total_predicted_lm_singlerun_constanttempfishing)$coefficients)]
#limit to interactions only (check this if there are any model changes!) row 2 and rows 34:64
model_coefs_reduced_predictions_singlerun_constanttempfishing <- model_coefs_reduced_predictions_singlerun_constanttempfishing[c(2,34:64),]
#adjust survey unit name by deleting beginning of string
model_coefs_reduced_predictions_singlerun_constanttempfishing[,survey_unit := substr(var, 17, str_length(var))][var == "year",survey_unit := "AI"]
#calculate interaction coefficients
AI_estimate <- model_coefs_reduced_predictions_singlerun_constanttempfishing[1,Estimate]
model_coefs_reduced_predictions_singlerun_constanttempfishing[1,estimate := AI_estimate]
model_coefs_reduced_predictions_singlerun_constanttempfishing[2:32,estimate := (AI_estimate + Estimate)]
predicted_dissim_trends_rnormruns_constanttempfishing <- rbind(predicted_dissim_trends_rnormruns_constanttempfishing, model_coefs_reduced_predictions_singlerun_constanttempfishing[,.(survey_unit, estimate)])
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90
[1] 91
[1] 92
[1] 93
[1] 94
[1] 95
[1] 96
[1] 97
[1] 98
[1] 99
[1] 100
[1] 101
[1] 102
[1] 103
[1] 104
[1] 105
[1] 106
[1] 107
[1] 108
[1] 109
[1] 110
[1] 111
[1] 112
[1] 113
[1] 114
[1] 115
[1] 116
[1] 117
[1] 118
[1] 119
[1] 120
[1] 121
[1] 122
[1] 123
[1] 124
[1] 125
[1] 126
[1] 127
[1] 128
[1] 129
[1] 130
[1] 131
[1] 132
[1] 133
[1] 134
[1] 135
[1] 136
[1] 137
[1] 138
[1] 139
[1] 140
[1] 141
[1] 142
[1] 143
[1] 144
[1] 145
[1] 146
[1] 147
[1] 148
[1] 149
[1] 150
[1] 151
[1] 152
[1] 153
[1] 154
[1] 155
[1] 156
[1] 157
[1] 158
[1] 159
[1] 160
[1] 161
[1] 162
[1] 163
[1] 164
[1] 165
[1] 166
[1] 167
[1] 168
[1] 169
[1] 170
[1] 171
[1] 172
[1] 173
[1] 174
[1] 175
[1] 176
[1] 177
[1] 178
[1] 179
[1] 180
[1] 181
[1] 182
[1] 183
[1] 184
[1] 185
[1] 186
[1] 187
[1] 188
[1] 189
[1] 190
[1] 191
[1] 192
[1] 193
[1] 194
[1] 195
[1] 196
[1] 197
[1] 198
[1] 199
[1] 200
[1] 201
[1] 202
[1] 203
[1] 204
[1] 205
[1] 206
[1] 207
[1] 208
[1] 209
[1] 210
[1] 211
[1] 212
[1] 213
[1] 214
[1] 215
[1] 216
[1] 217
[1] 218
[1] 219
[1] 220
[1] 221
[1] 222
[1] 223
[1] 224
[1] 225
[1] 226
[1] 227
[1] 228
[1] 229
[1] 230
[1] 231
[1] 232
[1] 233
[1] 234
[1] 235
[1] 236
[1] 237
[1] 238
[1] 239
[1] 240
[1] 241
[1] 242
[1] 243
[1] 244
[1] 245
[1] 246
[1] 247
[1] 248
[1] 249
[1] 250
[1] 251
[1] 252
[1] 253
[1] 254
[1] 255
[1] 256
[1] 257
[1] 258
[1] 259
[1] 260
[1] 261
[1] 262
[1] 263
[1] 264
[1] 265
[1] 266
[1] 267
[1] 268
[1] 269
[1] 270
[1] 271
[1] 272
[1] 273
[1] 274
[1] 275
[1] 276
[1] 277
[1] 278
[1] 279
[1] 280
[1] 281
[1] 282
[1] 283
[1] 284
[1] 285
[1] 286
[1] 287
[1] 288
[1] 289
[1] 290
[1] 291
[1] 292
[1] 293
[1] 294
[1] 295
[1] 296
[1] 297
[1] 298
[1] 299
[1] 300
[1] 301
[1] 302
[1] 303
[1] 304
[1] 305
[1] 306
[1] 307
[1] 308
[1] 309
[1] 310
[1] 311
[1] 312
[1] 313
[1] 314
[1] 315
[1] 316
[1] 317
[1] 318
[1] 319
[1] 320
[1] 321
[1] 322
[1] 323
[1] 324
[1] 325
[1] 326
[1] 327
[1] 328
[1] 329
[1] 330
[1] 331
[1] 332
[1] 333
[1] 334
[1] 335
[1] 336
[1] 337
[1] 338
[1] 339
[1] 340
[1] 341
[1] 342
[1] 343
[1] 344
[1] 345
[1] 346
[1] 347
[1] 348
[1] 349
[1] 350
[1] 351
[1] 352
[1] 353
[1] 354
[1] 355
[1] 356
[1] 357
[1] 358
[1] 359
[1] 360
[1] 361
[1] 362
[1] 363
[1] 364
[1] 365
[1] 366
[1] 367
[1] 368
[1] 369
[1] 370
[1] 371
[1] 372
[1] 373
[1] 374
[1] 375
[1] 376
[1] 377
[1] 378
[1] 379
[1] 380
[1] 381
[1] 382
[1] 383
[1] 384
[1] 385
[1] 386
[1] 387
[1] 388
[1] 389
[1] 390
[1] 391
[1] 392
[1] 393
[1] 394
[1] 395
[1] 396
[1] 397
[1] 398
[1] 399
[1] 400
[1] 401
[1] 402
[1] 403
[1] 404
[1] 405
[1] 406
[1] 407
[1] 408
[1] 409
[1] 410
[1] 411
[1] 412
[1] 413
[1] 414
[1] 415
[1] 416
[1] 417
[1] 418
[1] 419
[1] 420
[1] 421
[1] 422
[1] 423
[1] 424
[1] 425
[1] 426
[1] 427
[1] 428
[1] 429
[1] 430
[1] 431
[1] 432
[1] 433
[1] 434
[1] 435
[1] 436
[1] 437
[1] 438
[1] 439
[1] 440
[1] 441
[1] 442
[1] 443
[1] 444
[1] 445
[1] 446
[1] 447
[1] 448
[1] 449
[1] 450
[1] 451
[1] 452
[1] 453
[1] 454
[1] 455
[1] 456
[1] 457
[1] 458
[1] 459
[1] 460
[1] 461
[1] 462
[1] 463
[1] 464
[1] 465
[1] 466
[1] 467
[1] 468
[1] 469
[1] 470
[1] 471
[1] 472
[1] 473
[1] 474
[1] 475
[1] 476
[1] 477
[1] 478
[1] 479
[1] 480
[1] 481
[1] 482
[1] 483
[1] 484
[1] 485
[1] 486
[1] 487
[1] 488
[1] 489
[1] 490
[1] 491
[1] 492
[1] 493
[1] 494
[1] 495
[1] 496
[1] 497
[1] 498
[1] 499
[1] 500
[1] 501
[1] 502
[1] 503
[1] 504
[1] 505
[1] 506
[1] 507
[1] 508
[1] 509
[1] 510
[1] 511
[1] 512
[1] 513
[1] 514
[1] 515
[1] 516
[1] 517
[1] 518
[1] 519
[1] 520
[1] 521
[1] 522
[1] 523
[1] 524
[1] 525
[1] 526
[1] 527
[1] 528
[1] 529
[1] 530
[1] 531
[1] 532
[1] 533
[1] 534
[1] 535
[1] 536
[1] 537
[1] 538
[1] 539
[1] 540
[1] 541
[1] 542
[1] 543
[1] 544
[1] 545
[1] 546
[1] 547
[1] 548
[1] 549
[1] 550
[1] 551
[1] 552
[1] 553
[1] 554
[1] 555
[1] 556
[1] 557
[1] 558
[1] 559
[1] 560
[1] 561
[1] 562
[1] 563
[1] 564
[1] 565
[1] 566
[1] 567
[1] 568
[1] 569
[1] 570
[1] 571
[1] 572
[1] 573
[1] 574
[1] 575
[1] 576
[1] 577
[1] 578
[1] 579
[1] 580
[1] 581
[1] 582
[1] 583
[1] 584
[1] 585
[1] 586
[1] 587
[1] 588
[1] 589
[1] 590
[1] 591
[1] 592
[1] 593
[1] 594
[1] 595
[1] 596
[1] 597
[1] 598
[1] 599
[1] 600
[1] 601
[1] 602
[1] 603
[1] 604
[1] 605
[1] 606
[1] 607
[1] 608
[1] 609
[1] 610
[1] 611
[1] 612
[1] 613
[1] 614
[1] 615
[1] 616
[1] 617
[1] 618
[1] 619
[1] 620
[1] 621
[1] 622
[1] 623
[1] 624
[1] 625
[1] 626
[1] 627
[1] 628
[1] 629
[1] 630
[1] 631
[1] 632
[1] 633
[1] 634
[1] 635
[1] 636
[1] 637
[1] 638
[1] 639
[1] 640
[1] 641
[1] 642
[1] 643
[1] 644
[1] 645
[1] 646
[1] 647
[1] 648
[1] 649
[1] 650
[1] 651
[1] 652
[1] 653
[1] 654
[1] 655
[1] 656
[1] 657
[1] 658
[1] 659
[1] 660
[1] 661
[1] 662
[1] 663
[1] 664
[1] 665
[1] 666
[1] 667
[1] 668
[1] 669
[1] 670
[1] 671
[1] 672
[1] 673
[1] 674
[1] 675
[1] 676
[1] 677
[1] 678
[1] 679
[1] 680
[1] 681
[1] 682
[1] 683
[1] 684
[1] 685
[1] 686
[1] 687
[1] 688
[1] 689
[1] 690
[1] 691
[1] 692
[1] 693
[1] 694
[1] 695
[1] 696
[1] 697
[1] 698
[1] 699
[1] 700
[1] 701
[1] 702
[1] 703
[1] 704
[1] 705
[1] 706
[1] 707
[1] 708
[1] 709
[1] 710
[1] 711
[1] 712
[1] 713
[1] 714
[1] 715
[1] 716
[1] 717
[1] 718
[1] 719
[1] 720
[1] 721
[1] 722
[1] 723
[1] 724
[1] 725
[1] 726
[1] 727
[1] 728
[1] 729
[1] 730
[1] 731
[1] 732
[1] 733
[1] 734
[1] 735
[1] 736
[1] 737
[1] 738
[1] 739
[1] 740
[1] 741
[1] 742
[1] 743
[1] 744
[1] 745
[1] 746
[1] 747
[1] 748
[1] 749
[1] 750
[1] 751
[1] 752
[1] 753
[1] 754
[1] 755
[1] 756
[1] 757
[1] 758
[1] 759
[1] 760
[1] 761
[1] 762
[1] 763
[1] 764
[1] 765
[1] 766
[1] 767
[1] 768
[1] 769
[1] 770
[1] 771
[1] 772
[1] 773
[1] 774
[1] 775
[1] 776
[1] 777
[1] 778
[1] 779
[1] 780
[1] 781
[1] 782
[1] 783
[1] 784
[1] 785
[1] 786
[1] 787
[1] 788
[1] 789
[1] 790
[1] 791
[1] 792
[1] 793
[1] 794
[1] 795
[1] 796
[1] 797
[1] 798
[1] 799
[1] 800
[1] 801
[1] 802
[1] 803
[1] 804
[1] 805
[1] 806
[1] 807
[1] 808
[1] 809
[1] 810
[1] 811
[1] 812
[1] 813
[1] 814
[1] 815
[1] 816
[1] 817
[1] 818
[1] 819
[1] 820
[1] 821
[1] 822
[1] 823
[1] 824
[1] 825
[1] 826
[1] 827
[1] 828
[1] 829
[1] 830
[1] 831
[1] 832
[1] 833
[1] 834
[1] 835
[1] 836
[1] 837
[1] 838
[1] 839
[1] 840
[1] 841
[1] 842
[1] 843
[1] 844
[1] 845
[1] 846
[1] 847
[1] 848
[1] 849
[1] 850
[1] 851
[1] 852
[1] 853
[1] 854
[1] 855
[1] 856
[1] 857
[1] 858
[1] 859
[1] 860
[1] 861
[1] 862
[1] 863
[1] 864
[1] 865
[1] 866
[1] 867
[1] 868
[1] 869
[1] 870
[1] 871
[1] 872
[1] 873
[1] 874
[1] 875
[1] 876
[1] 877
[1] 878
[1] 879
[1] 880
[1] 881
[1] 882
[1] 883
[1] 884
[1] 885
[1] 886
[1] 887
[1] 888
[1] 889
[1] 890
[1] 891
[1] 892
[1] 893
[1] 894
[1] 895
[1] 896
[1] 897
[1] 898
[1] 899
[1] 900
[1] 901
[1] 902
[1] 903
[1] 904
[1] 905
[1] 906
[1] 907
[1] 908
[1] 909
[1] 910
[1] 911
[1] 912
[1] 913
[1] 914
[1] 915
[1] 916
[1] 917
[1] 918
[1] 919
[1] 920
[1] 921
[1] 922
[1] 923
[1] 924
[1] 925
[1] 926
[1] 927
[1] 928
[1] 929
[1] 930
[1] 931
[1] 932
[1] 933
[1] 934
[1] 935
[1] 936
[1] 937
[1] 938
[1] 939
[1] 940
[1] 941
[1] 942
[1] 943
[1] 944
[1] 945
[1] 946
[1] 947
[1] 948
[1] 949
[1] 950
[1] 951
[1] 952
[1] 953
[1] 954
[1] 955
[1] 956
[1] 957
[1] 958
[1] 959
[1] 960
[1] 961
[1] 962
[1] 963
[1] 964
[1] 965
[1] 966
[1] 967
[1] 968
[1] 969
[1] 970
[1] 971
[1] 972
[1] 973
[1] 974
[1] 975
[1] 976
[1] 977
[1] 978
[1] 979
[1] 980
[1] 981
[1] 982
[1] 983
[1] 984
[1] 985
[1] 986
[1] 987
[1] 988
[1] 989
[1] 990
[1] 991
[1] 992
[1] 993
[1] 994
[1] 995
[1] 996
[1] 997
[1] 998
[1] 999
[1] 1000
#reduce to mean and standard deviation
predicted_dissim_trends_rnormruns_constanttempfishing[,mean_dissim_coef:= mean(estimate),survey_unit][,sd_dissim := sd(estimate),.(survey_unit)]
predicted_dissim_trends_rnormrunsconstant_tempfishing.summary <- unique(predicted_dissim_trends_rnormruns_constanttempfishing[,.(survey_unit, mean_dissim_coef, sd_dissim)])
predicted_dissim_trends_rnormrunsconstant_tempfishing.summary[,pred_type := "fishing_and_temp_constant"]
predicted_dissim_trends_rnormruns.summary <- rbind(predicted_dissim_trends_rnormruns.summary, predicted_dissim_trends_rnormrunsconstant_tempfishing.summary)
Plotting observed vs predicted
jaccard_fishing_temp_model_observed_predicted_dt <- jaccard_total_coefs.r[predicted_dissim_trends_rnormruns.summary, on = "survey_unit"]
jaccard_fishing_temp_model_observed_predicted_dt[,pred_lower := mean_dissim_coef-sd_dissim][,pred_upper := mean_dissim_coef+sd_dissim]
#FULL MODEL, both temperature and fishing are allowed to vary
jaccard_fishing_temp_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"])
summary(jaccard_fishing_temp_model_observed_predicted_lm) #R^2 0.56
Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type ==
"full"])
Residuals:
Min 1Q Median 3Q Max
-0.0022236 -0.0006394 -0.0001630 0.0003886 0.0039875
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.548e-05 2.631e-04 0.097 0.923498
mean_dissim_coef 1.145e+00 3.086e-01 3.709 0.000843 ***
---
Signif. codes: 0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.001482 on 30 degrees of freedom
Multiple R-squared: 0.3144, Adjusted R-squared: 0.2916
F-statistic: 13.76 on 1 and 30 DF, p-value: 0.000843
(jaccard_fishing_temp_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "full"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
geom_point(aes(y = estimate, x = mean_dissim_coef)) +
geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
geom_abline(aes(slope = 1, intercept = 0)) +
lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n") +
theme_classic()
)
#fishing constant (fishing constant; temperature varies only)
jaccard_fishing_constant_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"])
summary(jaccard_fishing_constant_model_observed_predicted_lm)
Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type ==
"fishing_constant"])
Residuals:
Min 1Q Median 3Q Max
-2.810e-03 -1.001e-03 -2.156e-05 7.396e-04 2.963e-03
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.812e-06 2.656e-04 0.018 0.98567
mean_dissim_coef 2.149e+00 5.926e-01 3.626 0.00105 **
---
Signif. codes: 0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.001492 on 30 degrees of freedom
Multiple R-squared: 0.3047, Adjusted R-squared: 0.2816
F-statistic: 13.15 on 1 and 30 DF, p-value: 0.001054
#Temperature as a predictor, not fishing = R^2 = 0.40 (drop in 16% of variance explained when you lose fishing as predictor)
(jaccard_fishing_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
geom_point(aes(y = estimate, x = mean_dissim_coef)) +
geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
geom_abline(aes(slope = 1, intercept = 0)) +
lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(temperature varies fishing constant)") +
theme_classic()
)
#temp constant (fishing only; temperature constant)
jaccard_temperature_constant_model_observed_predicted_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"])
summary(jaccard_temperature_constant_model_observed_predicted_lm) #0.26 Drop in 30% of variance explained when you lose temperature as a predictor
Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type ==
"temp_constant"])
Residuals:
Min 1Q Median 3Q Max
-0.0023482 -0.0005425 -0.0002359 0.0003885 0.0042139
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.524e-05 2.743e-04 0.201 0.84174
mean_dissim_coef 9.948e-01 3.128e-01 3.180 0.00341 **
---
Signif. codes: 0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.001548 on 30 degrees of freedom
Multiple R-squared: 0.2521, Adjusted R-squared: 0.2272
F-statistic: 10.11 on 1 and 30 DF, p-value: 0.003408
(jaccard_temperature_constant_model_observed_predicted <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "temp_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
geom_point(aes(y = estimate, x = mean_dissim_coef)) +
geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm") +
geom_abline(aes(slope = 1, intercept = 0)) +
lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing varies temperature constant)") +
theme_classic()
)
#both temperature and fish held constant
jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm <- lm(estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"])
summary(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey_lm) #%20 #drop of 36 from full
Call:
lm(formula = estimate ~ mean_dissim_coef, data = jaccard_fishing_temp_model_observed_predicted_dt[pred_type ==
"fishing_and_temp_constant"])
Residuals:
Min 1Q Median 3Q Max
-0.0028055 -0.0008924 -0.0002188 0.0007432 0.0035791
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.548e-05 2.829e-04 0.196 0.84587
mean_dissim_coef 1.722e+00 6.186e-01 2.783 0.00922 **
---
Signif. codes: 0 ā***ā 0.001 ā**ā 0.01 ā*ā 0.05 ā.ā 0.1 ā ā 1
Residual standard error: 0.001595 on 30 degrees of freedom
Multiple R-squared: 0.2052, Adjusted R-squared: 0.1787
F-statistic: 7.747 on 1 and 30 DF, p-value: 0.009223
(jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey <- ggplot(jaccard_fishing_temp_model_observed_predicted_dt[pred_type == "fishing_and_temp_constant"]) +
geom_errorbar(aes(x = mean_dissim_coef, ymin = lwr, ymax = upr), color = "lightgrey", linewidth = 0.4) +
geom_errorbarh(aes(y = estimate, xmin = mean_dissim_coef-sd_dissim, xmax = mean_dissim_coef+sd_dissim), color = "lightgrey", linewidth = 0.4) +
geom_point(aes(y = estimate, x = mean_dissim_coef)) +
geom_smooth(aes(y = estimate, x = mean_dissim_coef), color = "darkgrey",linetype = "dotted", method = "lm")+
geom_abline(aes(slope = 1, intercept = 0)) +
lims(x = c(min(jaccard_fishing_temp_model_observed_predicted_dt$pred_lower),max(jaccard_fishing_temp_model_observed_predicted_dt$pred_upper))) +
labs(y = "Observed β-diversity trend",x = "Predicted β-diversity trend\n(fishing and temperature constant)") +
theme_classic()
)
#merge
jaccard_fishing_sbt_model_observed_predicted_merge <- plot_grid(jaccard_fishing_temp_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
jaccard_fishing_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
jaccard_temperature_constant_model_observed_predicted + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")),
jaccard_fishing_temp_model_observed_predicted_tempfishconstantinsurvey + theme(plot.margin = unit(c(0.1,0.3,0.1,0.1),"cm")), ncol = 2, labels = c("a.","b.","c.","d."))
`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'`geom_smooth()` using formula = 'y ~ x'
ggsave(jaccard_fishing_sbt_model_observed_predicted_merge, path = here::here("figures"),filename = "jaccard_fishing_sbt_model_observed_predicted_merge.jpg", height =6, width = 8)
NA
NA
###Letās visualize model coefficients with temperature and fishing (similar to figure 2)
#extract coefficients using linear algebra
interaction_avg_model_coef <- data.table(lm_interaction_coefficients_se(mod_name = model_avg_delta4, model_avg = T, SBT_fish = T))
Warning: NaNs produced
#add survey names
interaction_avg_model_coef[,survey_unit := factor(c(all_surveys[!all_surveys %in% c("ROCKALL", "GSL-S")],all_surveys[!all_surveys %in% c("ROCKALL", "GSL-S")]))]
#add predictor names
interaction_avg_model_coef[,predictor := c(rep("Relative fishing catch",32),rep("Minimum temperature",32))]
#reorder temperature and fishing
interaction_avg_model_coef[,predictor := factor(predictor, levels = c("Minimum temperature","Relative fishing catch"))]
#link for full survey name
interaction_avg_model_coef <- color_link[interaction_avg_model_coef, on = "survey_unit"]
#list of Survey Name Season in order by survey unit
setkey(interaction_avg_model_coef, survey_unit)
survey_name_season_ordered <- unique(interaction_avg_model_coef$Survey_Name_Season)
#Make Survey_Name_Season a factor
interaction_avg_model_coef[,Survey_Name_Season := factor(Survey_Name_Season, levels = survey_name_season_ordered)]
#mark significance
interaction_avg_model_coef[,Significant := ifelse((estimate-se > 0 & estimate+se > 0) | (estimate-se < 0 & estimate+se < 0),T,F)]
#Plot both
sbt_fishing_avg_model_coef <- ggplot() +
geom_hline(yintercept = 0, linetype = "dotted") +
geom_point(data = interaction_avg_model_coef, aes(x = Survey_Name_Season, y = estimate, color = Significant)) +
geom_errorbar(data = interaction_avg_model_coef, aes(x = Survey_Name_Season, ymin = estimate-se, ymax = estimate+se, color = Significant), width = 0) +
scale_color_manual(values = c("darkgrey","black")) +
facet_wrap(~predictor, scales = "free_x") +
scale_x_discrete(limits = rev) +
labs(y = "Coefficient", x = "") +
coord_flip() +
theme_classic()
Plot all other coefficients in averaged model
model_avg_values.nonfishortemp <- model_avg_values[c(2,6,102:106),]
#mark significance
model_avg_values.nonfishortemp[,Significant := ifelse((coef-se > 0 & coef+se > 0) | (coef-se < 0 & coef+se < 0),T,F)]
#more helpful names for variables
model_avg_values.nonfishortemp[,Variable := c("Area","Species count","Depth","Number of tows","Latitude","Depth range","Latitude range")]
#make factor with order
model_avg_values.nonfishortemp[,Variable := factor(Variable, levels = c("Area","Species count","Number of tows","Depth","Depth range","Latitude","Latitude range"))]
#plot
all_avg_model_coef <- ggplot() +
geom_point(data = model_avg_values.nonfishortemp, aes(x = Variable, y = coef, color = Significant)) +
geom_errorbar(data = model_avg_values.nonfishortemp, aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
labs(y = "Coefficient", x = "\n\n\n") +
coord_flip() +
theme_classic()
#plot
all_but_lat_avg_model_coef <- ggplot() +
geom_point(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, y = coef, color = Significant)) +
geom_errorbar(data = model_avg_values.nonfishortemp[Variable == "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()
#plot
lat_avg_model_coef <- ggplot() +
geom_point(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, y = coef, color = Significant)) +
geom_errorbar(data = model_avg_values.nonfishortemp[Variable != "Latitude"], aes(x = Variable, ymin = coef-se, ymax = coef+se, color = Significant), width = 0) +
scale_color_manual(values = c("darkgrey","black")) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev) +
labs(y = "Coefficient", x = "Model variable") +
coord_flip() +
theme_classic()
#merge into single plot
model_coef_summary_sbt_jaccard <- cowplot::plot_grid(sbt_fishing_avg_model_coef+theme(legend.position = "null", axis.title.x = element_blank()),all_avg_model_coef+theme(legend.position = "null"), ncol = 1, labels = c(" a. b."," c."), label_y = 0.99, rel_heights = c(3,1))
#Figure 3
ggsave(model_coef_summary_sbt_jaccard, path = here::here("figures"),filename = "model_coef_summary_sbt_jaccard.jpg", height = 6.5, width = 8, unit = "in")
ggsave(model_coef_summary_sbt_jaccard, path = here::here("figures"),filename = "Fig3.tiff", height = 6.5, width = 8, unit = "in", dpi = 300, compression = "lzw")
###Investigating whatās going on with Aleutian Islands
#DISSIMILARITY VS TEMPERATURE
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = yearly_min_bypoint_avg, y = annual_dissimilarity_value)) +
theme_classic()
#If we ignore two coldest years, we have roughly points at (2.07,0.76) and (2.2, 0.69), slope = -1.8, which matches what we get (also note large confidence bounds)
#So, the question is, what is leading the model to ignore these two years? What explains this variation in the coldest years?
#DISSIMILARITY
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = area_km, y = annual_dissimilarity_value)) +
theme_classic()
#Very large sampling area and very low sampling areas are both lowest values
#TEMPERATURE
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = area_km, y = yearly_min_bypoint_avg)) +
theme_classic()
#Very large sampling area coincides with lowest temperatures
#DISSIMILARITY
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = spp_count_annual, y = annual_dissimilarity_value)) +
theme_classic()
#Annual dissimilarity value increases with richness
#TEMPERATURE
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = spp_count_annual, y = yearly_min_bypoint_avg)) +
theme_classic()
#Lowest temperature observation coincides with lowest richness
#DISSIMILARITY
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = haul_id_count_annual, y = annual_dissimilarity_value)) +
theme_classic()
#Low number of hauls led to very low annual dissimilarity values
#TEMPERATURE
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = haul_id_count_annual, y = yearly_min_bypoint_avg)) +
theme_classic()
#Lowest temperature matched with lowest number of samples
#DISSIMILARITY
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = depth_annual_avg, y = annual_dissimilarity_value)) +
theme_classic()
#Some deep trawls have very low annual dissimilarity, while other deep trawls have high annual dissimilarity
#TEMPERATURE
ggplot(dissimilarity_covariates_dredge.dt[survey_unit == "AI"]) +
geom_point(aes(x = depth_annual_avg, y = yearly_min_bypoint_avg)) +
theme_classic()
#Deepest trawls happened on cold days
What leads to variance being explained by variables OTHER than temperature in coldest years So, in summary - Very large sampling area and very small sampling areas are both linked to years of low dissimilarity - Very large sampling area coincides with lowest temperatures - Annual dissimilarity value increases with richness - Lowest temperature observation coincides with lowest richness - Low number of hauls led to very low annual dissimilarity values - Lowest temperature matched with lowest number of samples - Some deep trawls have very low annual dissimilarity, while other deep trawls have high annual dissimilarity - Deepest trawls happened on cold days
For the Aleutian islands, the coldest years were also: - the years of largest sampling area - the years of lowest richness - the years of lowest sampling effort (# tows) - the years of deepest sampling
#IF WE ACTUALLY USE THIS, NEED TO UPDATE ##Show intercepts by region
and season #
#{r} # ##survey intercepts #survey_intercepts <- model_avg_values[c(1,8:38),.(coef_name, coef, se)] #survey_intercepts[1,survey_unit := "AI"][1,coef_true := coef] #survey_intercepts[2:32,survey_unit := substr(coef_name, 12, str_length(coef_name))][2:32,coef_true := survey_intercepts[1,coef_true]+coef] #survey_intercepts <- color_link[survey_intercepts, on = "survey_unit"] #survey_intercepts[,Survey_Name_Season := reorder(Survey_Name_Season, coef_true)] # # ##season_intercepts #season_intercepts <- model_avg_values[c(1,3:5),.(coef_name, coef, se)] #season_intercepts[1,season := "Spring"][1,coef_true := coef] #season_intercepts[2:4,season := substr(coef_name, 7, str_length(coef_name))][2:4,coef_true := season_intercepts[1,coef_true]+coef] #season_intercepts[,season := reorder(season, coef_true)] # # ##survey intercepts #survey_model_coef <- ggplot() + #geom_point(data = survey_intercepts, aes(x = Survey_Name_Season, y = coef_true)) + #geom_errorbar(data = survey_intercepts, aes(x = Survey_Name_Season, ymin = coef_true-se, ymax = coef_true+se), width = 0) + #scale_x_discrete(limits = rev) + ## ylim(0.35,1) + # labs(y = "Intercept", x = "") + #coord_flip() + #theme_classic() # ##season intercepts #season_model_coef <- ggplot() + #geom_point(data = season_intercepts, aes(x = season, y = coef_true)) + #geom_errorbar(data = season_intercepts, aes(x = season, ymin = coef_true-se, ymax = coef_true+se), width = 0) + #scale_x_discrete(limits = rev) + # # ylim(0.35,1) + # labs(y = "Intercept", x = "") + #coord_flip() + #theme_classic() # ##merge into single plot # #model_intercept_jaccard <- cowplot::plot_grid(survey_model_coef+theme(axis.title.x = element_blank()), # season_model_coef, ncol = 1, labels = c("a.","b."), label_y = 0.99, rel_heights = c(3,1), align = "v") # #ggsave(model_intercept_jaccard, path = here::here("figures"),filename = "model_intercept_jaccard.jpg", height = 6.5, width = 6.5, unit = "in") # # # ##